! This is a 1D advection example using square initial condition and periodic
! boundary condition for TSPAS finite difference scheme.
!
! Li Dong <dongli@lasg.iap.ac.cn>
!
! - 2018-03-21: Initial creation.
! - 2024-03-24: Use adv_1d_test_case_mod.

program tspas_adv_1d_case

  use adv_1d_test_case_mod

  implicit none

  real, allocatable :: rho  (:,:)   ! Tracer density being advected at cell centers
  real, allocatable :: flx  (:)     ! Flux at cell interfaces
  real, allocatable :: gamma(:)     ! cfl * (1 - cfl) at cell interfaces
  real, allocatable :: A    (:)     ! 
  real, allocatable :: cstar(:)     ! Switch of upwind scheme
  real, allocatable :: ustar(:)     ! Modified velocity
  real coef                         ! dt / dx
  real, parameter :: eps = 1.0d-80  ! A small value to avoid divided-by-zero
  integer, parameter :: ns = 1      ! Stencil width
  integer, parameter :: star = 0
  integer i
  character(30), parameter :: scheme = 'tspas'

  namelist /params/ nx, nt, dt, u

  call get_command_argument(1, namelist_path)
  inquire(file=namelist_path, exist=is_exist)
  if (is_exist) then
    open(10, file=namelist_path)
    read(10, nml=params)
    close(10)
  end if

  allocate(rho  (1-ns:nx+ns,0:2))
  allocate(flx  (1-ns:nx+ns))
  allocate(gamma(1-ns:nx+ns))
  allocate(A    (1-ns:nx+ns))
  allocate(cstar(1-ns:nx+ns))
  allocate(ustar(1-ns:nx+ns))

  call adv_1d_test_case_init('square', ns, rho(:,old))
  call output(scheme, 0, ns, nx, x, rho(:,old))

  ! Run integration.
  print *, time_step, sum(rho(1:nx,old))
  do while (time_step < nt)
    call tspas(rho(:,old), rho(:,star), flx)
    do i = 1, nx
      rho(i,new) = rho(i,old) - dt / dx * (flx(i) - flx(i-1))
    end do
    call apply_bc(ns, nx, rho(:,new))
    call advance_time()
    call output(scheme, time_step, ns, nx, x, rho(:,old))
    print *, time_step, sum(rho(1:nx,old))
  end do

  deallocate(rho, flx, gamma, A, cstar, ustar)

  call adv_1d_test_case_final()

contains

  subroutine tspas(q, qstar, f)

    real, intent(in   ) :: q    (1-ns:nx+ns)
    real, intent(inout) :: qstar(1-ns:nx+ns)
    real, intent(out  ) :: f    (1-ns:nx+ns)

    real c, alpha, beta, tmp1, tmp2
    integer i

    c = dt / dx
    ! Run Lax-Wendroff pass.
    do i = 1, nx
      f(i) = 0.5d0 * (u * (q(i+1) + q(i)) - c * u**2 * (q(i+1) - q(i)))
    end do
    call apply_bc(ns, nx, f)
    do i = 1, nx + 1
      alpha = abs(u) * coef
      gamma(i) = alpha * (1.0d0 - alpha)
    end do
    do i = 1, nx
      beta = max(2.0d0 / (2.0d0 - gamma(i)), 2.0d0 / (2.0d0 - gamma(i+1)))
      qstar(i) = q(i) - beta * c * (f(i) - f(i-1))
    end do
    call apply_bc(ns, nx, qstar)
    ! Calculate A.
    do i = 1, nx
      A(i) = (qstar(i) - max(q(i-1), q(i), q(i+1))) * &
             (qstar(i) - min(q(i-1), q(i), q(i+1)))
    end do
    call apply_bc(ns, nx, A)
    ! Calculate u_star
    do i = 1, nx
      tmp1 = abs(A(i  )) + A(i  )
      tmp2 = abs(A(i+1)) + A(i+1)
      cstar(i) = 0.5d0  * (tmp1 / (abs(A(i)) + eps) + tmp2 / (abs(A(i+1)) + eps)) - &
                 0.25d0 * (tmp1 * tmp2 / (abs(A(i) * A(i+1)) + eps))
      ustar(i) = (cstar(i) + (1 - cstar(i)) * dt / dx * abs(u)) * u
    end do
    ! Run upwind pass.
    do i = 1, nx
      f(i) = 0.5d0 * (u * (q(i+1) + q(i)) - abs(ustar(i)) * (q(i+1) - q(i)))
    end do
    call apply_bc(ns, nx, f)

  end subroutine tspas

end program tspas_adv_1d_case
